Load all required libraries.

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.3     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(broom)

Read in raw data from RDS.

raw_data <- readRDS("./year2.RDS")

Make a few small modifications to names and data for visualizations.

final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
  rename(Facility = wrf) %>%
  mutate(Facility = recode(Facility, 
                           "NO" = "WRF A",
                           "MI" = "WRF B",
                           "CC" = "WRF C"))

Seperate the data by gene target to ease layering in the final plot

#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>% 
  select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke)) %>%
  group_by(date) %>% summarise_if(is.numeric, mean)

#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')
#remove facilty C for now
#only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
#only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]

only_n1 <- only_n1[!(only_n1$Facility == "WRF A" & only_n1$date == "2020-11-02"), ]
only_n2 <- only_n2[!(only_n2$Facility == "WRF A" & only_n2$date == "2020-11-02"), ]

Build the main plot

      #first layer is the background epidemic curve
        p1 <- only_background %>%
              plotly::plot_ly() %>%
              plotly::add_trace(x = ~date, y = ~new_cases_clarke, 
                                type = "bar", 
                                hoverinfo = "text",
                                text = ~paste('</br> Date: ', date,
                                                     '</br> Daily Cases: ', new_cases_clarke),
                                alpha = 0.5,
                                name = "Daily Reported Cases",
                                color = background_color,
                                colors = background_color,
                                showlegend = FALSE) %>%
            layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #renders the main plot layer two as seven day moving average
        p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke, 
                             type = "scatter",
                             mode = "lines",
                             hoverinfo = "text",
                            text = ~paste('</br> Date: ', date,
                                                     '</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
                             name = "Seven Day Moving Average Athens",
                             line = list(color = seven_day_ave_color),
                             showlegend = FALSE)
      

        
        #renders the main plot layer three as positive target hits
        
        p2 <- plotly::plot_ly() %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n1,
                                       symbol = ~Facility,
                                       marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n2,
                                       symbol = ~Facility,
                                       marker = list(color = '#D95F02', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
            layout(yaxis = list(title = "SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 type = "log",
                                 dtick = 1,
                                 automargin = TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #adds the limit of detection dashed line
        p2 <- p2 %>% plotly::add_segments(x = as.Date("2021-06-30"), 
                                          xend = ~max(date + 10), 
                                          y = 3571.429, yend = 3571.429,
                                          opacity = 0.35,
                                          line = list(color = "black", dash = "dash")) %>%
          layout(annotations = list(x = as.Date("2021-06-30"), y = 3.8, xref = "x", yref = "y", 
                                    text = "Limit of Detection", showarrow = FALSE))

        

        p1
        p2

Combine the two main plot pieces as a subplot

#seperate n1 and n2 frames by site
#n1
wrf_a_only_n1 <- subset(only_n1, Facility == "WRF A")
wrf_b_only_n1 <- subset(only_n1, Facility == "WRF B")
wrf_c_only_n1 <- subset(only_n1, Facility == "WRF C")

#n2
wrf_a_only_n2 <- subset(only_n2, Facility == "WRF A")
wrf_b_only_n2 <- subset(only_n2, Facility == "WRF B")
wrf_c_only_n2 <- subset(only_n2, Facility == "WRF C")


#rejoin the old data frames then seperate in to averages for each plant. 
wrfa_both <- full_join(wrf_a_only_n1, wrf_a_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
wrfb_both <- full_join(wrf_b_only_n1, wrf_b_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
wrfc_both <- full_join(wrf_c_only_n1, wrf_c_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
#get max date
maxdate <- max(wrfa_both$date)
mindate <- min(wrfa_both$date)

Build loess smoothing figures figures

This makes the individual plots

#**************************************WRF A PLOT**********************************************
#add trendlines 
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_botha <- ggplot(wrfa_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_botha<<-..y..), method = "loess", color = '#1B9E77', 
              span = 0.3, n = 260)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_botha
## `geom_smooth()` using formula 'y ~ x'

fit_botha
##   [1] 11.53127 11.56979 11.60796 11.64580 11.68331 11.72049 11.75734 11.79389
##   [9] 11.83012 11.86605 11.90168 11.93701 11.97205 12.00680 12.04122 12.07533
##  [17] 12.10912 12.14260 12.17575 12.20858 12.24123 12.27377 12.30618 12.33842
##  [25] 12.37043 12.40217 12.43361 12.46470 12.49539 12.52566 12.55544 12.58470
##  [33] 12.61341 12.64150 12.66975 12.69869 12.72796 12.75721 12.78607 12.81419
##  [41] 12.84121 12.86678 12.89053 12.91502 12.94209 12.97028 12.99812 13.02415
##  [49] 13.04691 13.06495 13.08024 13.09556 13.11061 13.12510 13.13876 13.15130
##  [57] 13.16243 13.17186 13.17931 13.18450 13.18713 13.18693 13.18360 13.17686
##  [65] 13.16544 13.14876 13.12758 13.10265 13.07472 13.04454 13.01288 12.98048
##  [73] 12.94811 12.91651 12.88643 12.85864 12.82810 12.79097 12.74983 12.70724
##  [81] 12.66576 12.62797 12.59641 12.56830 12.53934 12.50982 12.48003 12.45025
##  [89] 12.42078 12.39191 12.36392 12.33712 12.31179 12.28822 12.26670 12.24753
##  [97] 12.23098 12.21772 12.20775 12.20039 12.19497 12.19082 12.18726 12.18362
## [105] 12.17922 12.17339 12.16744 12.16294 12.15960 12.15712 12.15521 12.15359
## [113] 12.15195 12.15059 12.14997 12.15003 12.15073 12.15204 12.15391 12.15629
## [121] 12.15914 12.16242 12.16609 12.17010 12.17440 12.17896 12.18374 12.18965
## [129] 12.19740 12.20659 12.21682 12.22772 12.23887 12.24988 12.26037 12.26994
## [137] 12.27820 12.28474 12.28918 12.29206 12.29430 12.29615 12.29783 12.29958
## [145] 12.30162 12.30419 12.30577 12.30499 12.30231 12.29817 12.29303 12.28733
## [153] 12.28155 12.27611 12.27149 12.26812 12.26647 12.26698 12.27011 12.27631
## [161] 12.28603 12.29973 12.31793 12.34039 12.36637 12.39514 12.42596 12.45810
## [169] 12.49083 12.52341 12.55510 12.58518 12.61291 12.63755 12.66471 12.69949
## [177] 12.74044 12.78612 12.83505 12.88580 12.93690 12.98690 13.03434 13.07778
## [185] 13.11576 13.14682 13.16951 13.18237 13.18815 13.19066 13.19003 13.18641
## [193] 13.17994 13.17075 13.15899 13.14480 13.12832 13.10968 13.08904 13.06652
## [201] 13.04227 13.01643 12.98408 12.94183 12.89211 12.83735 12.77999 12.72247
## [209] 12.66721 12.61666 12.57324 12.53119 12.48436 12.43462 12.38385 12.33393
## [217] 12.28673 12.24414 12.20403 12.16319 12.12187 12.08030 12.03873 11.99739
## [225] 11.95653 11.91638 11.87719 11.83919 11.80263 11.76775 11.73478 11.70397
## [233] 11.67459 11.64585 11.61785 11.59072 11.56458 11.53954 11.51572 11.49324
## [241] 11.47222 11.45278 11.43503 11.41909 11.40509 11.39278 11.38186 11.37232
## [249] 11.36417 11.35741 11.35204 11.34805 11.34546 11.34426 11.34444 11.34602
## [257] 11.34899 11.35335 11.35911 11.36626
#assign fits to a vector
both_trenda <- fit_botha

#extract y min and max for each
limits_botha <- ggplot_build(extract_botha)$data
## `geom_smooth()` using formula 'y ~ x'
limits_botha <- as.data.frame(limits_botha)
both_ymina <- limits_botha$ymin
both_ymaxa <- limits_botha$ymax

#reassign dataframes (just to be safe)
work_botha <- wrfa_both

#fill in missing dates to smooth fits
work_botha <- work_botha %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_botha <- work_botha$date

#create a new smooth dataframe to layer
smooth_frame_botha <- data.frame(date_vec_botha, both_trenda, both_ymina, both_ymaxa)
#WRF A
#plot smooth frames
p_wrf_a <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_botha, y = ~both_trenda,
                    data = smooth_frame_botha,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_botha,
                                  '</br> Median Log Copies: ', round(both_trenda, digits = 2)),
                    line = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_botha, ymin = ~both_ymina, ymax = ~both_ymaxa,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_botha, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxa, digits = 2),
                                  '</br> Min Log Copies: ', round(both_ymina, digits = 2)),
                    name = "",
                    fillcolor = '#1B9E77',
                    line = list(color = '#1B9E77')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF A") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfa_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#1B9E77', size = 6, opacity = 0.65))

p_wrf_a
save(p_wrf_a, file = "./site_objects/wrf_a_year2.rda")
#**************************************WRF B PLOT**********************************************
#add trendlines 
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothb <- ggplot(wrfb_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_bothb<<-..y..), method = "loess", color = '#D95F02', 
              span = 0.3, n = 260)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_bothb
## `geom_smooth()` using formula 'y ~ x'

fit_bothb
##   [1] 10.71297 10.80589 10.89692 10.98605 11.07327 11.15857 11.24192 11.32332
##   [9] 11.40275 11.48020 11.55565 11.62910 11.70052 11.76990 11.83726 11.90263
##  [17] 11.96605 12.02756 12.08721 12.14503 12.20079 12.25428 12.30556 12.35470
##  [25] 12.40178 12.44687 12.49003 12.53134 12.57087 12.60868 12.64485 12.67945
##  [33] 12.71255 12.74421 12.77384 12.80097 12.82581 12.84862 12.86963 12.88908
##  [41] 12.90722 12.92426 12.94047 12.95365 12.96231 12.96767 12.97096 12.97341
##  [49] 12.97625 12.98071 12.98535 12.98804 12.98903 12.98853 12.98678 12.98402
##  [57] 12.98047 12.97637 12.97194 12.96742 12.96303 12.95902 12.95561 12.95303
##  [65] 12.95097 12.94888 12.94666 12.94419 12.94136 12.93805 12.93415 12.92955
##  [73] 12.92414 12.91779 12.91041 12.90187 12.89328 12.88545 12.87781 12.86976
##  [81] 12.86075 12.85017 12.83747 12.82300 12.80757 12.79125 12.77412 12.75626
##  [89] 12.73774 12.71864 12.69904 12.67901 12.65864 12.63799 12.61714 12.59618
##  [97] 12.57517 12.55382 12.53174 12.50886 12.48515 12.46054 12.43499 12.40844
## [105] 12.38083 12.35212 12.31765 12.27508 12.22782 12.17928 12.13287 12.09200
## [113] 12.06009 12.03236 12.00226 11.97055 11.93797 11.90527 11.87320 11.84251
## [121] 11.81395 11.78826 11.76620 11.74851 11.73594 11.72924 11.72917 11.73598
## [129] 11.74880 11.76668 11.78866 11.81378 11.84108 11.86961 11.89841 11.92652
## [137] 11.95298 11.97683 11.99713 12.02054 12.05234 12.08947 12.12884 12.16738
## [145] 12.20203 12.22971 12.25357 12.27870 12.30491 12.33202 12.35984 12.38820
## [153] 12.41691 12.44578 12.47464 12.50330 12.53158 12.55929 12.58626 12.61229
## [161] 12.63721 12.66083 12.68322 12.70467 12.72529 12.74519 12.76450 12.78334
## [169] 12.80181 12.82005 12.83816 12.85627 12.87449 12.89295 12.91412 12.93983
## [177] 12.96921 13.00137 13.03544 13.07055 13.10582 13.14038 13.17334 13.20384
## [185] 13.23099 13.25393 13.27176 13.28363 13.29152 13.29796 13.30291 13.30629
## [193] 13.30806 13.30817 13.30655 13.30314 13.29790 13.29076 13.28168 13.27058
## [201] 13.25742 13.24215 13.22424 13.20346 13.18018 13.15474 13.12750 13.09882
## [209] 13.06905 13.03854 13.00767 12.97619 12.94359 12.90980 12.87473 12.83832
## [217] 12.80049 12.76115 12.71822 12.67025 12.61812 12.56272 12.50491 12.44560
## [225] 12.38565 12.32596 12.26740 12.21085 12.15720 12.10733 12.06212 12.02245
## [233] 11.98617 11.95055 11.91563 11.88144 11.84801 11.81538 11.78359 11.75266
## [241] 11.72263 11.69354 11.66543 11.63831 11.61224 11.58718 11.56312 11.54006
## [249] 11.51800 11.49695 11.47693 11.45792 11.43996 11.42303 11.40715 11.39232
## [257] 11.37856 11.36587 11.35425 11.34372
#assign fits to a vector
both_trendb <- fit_bothb

#extract y min and max for each
limits_bothb <- ggplot_build(extract_bothb)$data
## `geom_smooth()` using formula 'y ~ x'
limits_bothb <- as.data.frame(limits_bothb)
both_yminb <- limits_bothb$ymin
both_ymaxb <- limits_bothb$ymax

#reassign dataframes (just to be safe)
work_bothb <- wrfb_both

#fill in missing dates to smooth fits
work_bothb <- work_bothb %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothb <- work_bothb$date

#create a new smooth dataframe to layer
smooth_frame_bothb <- data.frame(date_vec_bothb, both_trendb, both_yminb, both_ymaxb)
#WRF B
#plot smooth frames
p_wrf_b <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_bothb, y = ~both_trendb,
                    data = smooth_frame_bothb,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothb,
                                  '</br> Median Log Copies: ', round(both_trendb, digits = 2)),
                    line = list(color = '#D95F02', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothb, ymin = ~both_yminb, ymax = ~both_ymaxb,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothb, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxb, digits = 2),
                                  '</br> Min Log Copies: ', round(both_yminb, digits = 2)),
                    name = "",
                    fillcolor = '#D95F02',
                    line = list(color = '#D95F02')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF B") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfb_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#D95F02', size = 6, opacity = 0.65))

p_wrf_b
save(p_wrf_b, file = "./site_objects/wrf_b_year2.rda")

#**************************************WRF C PLOT********************************************** #add trendlines #extract data from geom_smooth # *********************************span 0.6*********************************** #*****************Must always update the n = TOTAL NUMBER OF DAYS*************************

extract_bothc <- ggplot(wrfc_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_bothc<<-..y..), method = "loess", color = '#E7298A', 
              span = 0.3, n = 260)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_bothc
## `geom_smooth()` using formula 'y ~ x'

fit_bothc
##   [1] 10.82125 10.88662 10.95077 11.01366 11.07528 11.13558 11.19454 11.25214
##   [9] 11.30833 11.36310 11.41641 11.46823 11.51853 11.56734 11.61475 11.66080
##  [17] 11.70554 11.74903 11.79132 11.83246 11.87188 11.90910 11.94428 11.97760
##  [25] 12.00924 12.03936 12.06814 12.09575 12.12237 12.14817 12.17332 12.19799
##  [33] 12.22236 12.24660 12.26998 12.29182 12.31238 12.33194 12.35076 12.36911
##  [41] 12.38726 12.40548 12.42402 12.44434 12.46690 12.49056 12.51416 12.53655
##  [49] 12.55657 12.57307 12.58711 12.60046 12.61304 12.62476 12.63552 12.64525
##  [57] 12.65384 12.66121 12.66727 12.67193 12.67510 12.67670 12.67662 12.67479
##  [65] 12.67088 12.66479 12.65679 12.64713 12.63605 12.62383 12.61070 12.59692
##  [73] 12.58275 12.56844 12.55424 12.54041 12.52540 12.50794 12.48870 12.46835
##  [81] 12.44755 12.42698 12.40730 12.38684 12.36382 12.33873 12.31205 12.28427
##  [89] 12.25588 12.22737 12.19923 12.17194 12.14600 12.12190 12.10011 12.08114
##  [97] 12.06546 12.05236 12.04065 12.03015 12.02065 12.01197 12.00391 11.99628
## [105] 11.98888 11.98153 11.97624 11.97430 11.97447 11.97550 11.97613 11.97512
## [113] 11.97122 11.96583 11.96114 11.95709 11.95361 11.95062 11.94806 11.94586
## [121] 11.94395 11.94227 11.94073 11.93928 11.93784 11.93634 11.93472 11.93292
## [129] 11.93091 11.92869 11.92622 11.92351 11.92053 11.91727 11.91372 11.90985
## [137] 11.90566 11.90113 11.89624 11.88731 11.87238 11.85405 11.83491 11.81755
## [145] 11.80456 11.79854 11.79672 11.79478 11.79292 11.79134 11.79024 11.78984
## [153] 11.79033 11.79193 11.79484 11.79927 11.80541 11.81348 11.82368 11.83622
## [161] 11.85130 11.86912 11.89060 11.91607 11.94488 11.97638 12.00991 12.04481
## [169] 12.08045 12.11615 12.15127 12.18516 12.21716 12.24662 12.27877 12.31839
## [177] 12.36414 12.41469 12.46868 12.52478 12.58166 12.63797 12.69237 12.74352
## [185] 12.79009 12.83074 12.86411 12.88889 12.90807 12.92549 12.94108 12.95474
## [193] 12.96640 12.97596 12.98335 12.98849 12.99127 12.99164 12.98949 12.98474
## [201] 12.97731 12.96712 12.95259 12.93284 12.90885 12.88158 12.85199 12.82106
## [209] 12.78975 12.75902 12.72984 12.69866 12.66235 12.62246 12.58058 12.53829
## [217] 12.49715 12.45875 12.41988 12.37682 12.33036 12.28132 12.23049 12.17868
## [225] 12.12670 12.07534 12.02540 11.97770 11.93303 11.89220 11.85600 11.82525
## [233] 11.79800 11.77179 11.74666 11.72264 11.69977 11.67809 11.65762 11.63842
## [241] 11.62052 11.60394 11.58873 11.57493 11.56256 11.55160 11.54201 11.53377
## [249] 11.52689 11.52138 11.51723 11.51445 11.51303 11.51299 11.51431 11.51701
## [257] 11.52108 11.52653 11.53335 11.54156
#assign fits to a vector
both_trendc <- fit_bothc

#extract y min and max for each
limits_bothc <- ggplot_build(extract_bothc)$data
## `geom_smooth()` using formula 'y ~ x'
limits_bothc <- as.data.frame(limits_bothc)
both_yminc <- limits_bothc$ymin
both_ymaxc <- limits_bothc$ymax

#reassign dataframes (just to be safe)
work_bothc <- wrfc_both

#fill in missing dates to smooth fits
work_bothc <- work_bothc %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothc <- work_bothc$date

#create a new smooth dataframe to layer
smooth_frame_bothc <- data.frame(date_vec_bothc, both_trendc, both_yminc, both_ymaxc)
#WRF C
#plot smooth frames
p_wrf_c <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_bothc, y = ~both_trendc,
                    data = smooth_frame_bothc,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothc,
                                  '</br> Median Log Copies: ', round(both_trendc, digits = 2)),
                    line = list(color = '#E7298A', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothc, ymin = ~both_yminc, ymax = ~both_ymaxc,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothc, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxc, digits = 2),
                                  '</br> Min Log Copies: ', round(both_yminc, digits = 2)),
                    name = "",
                    fillcolor = '#E7298A',
                    line = list(color = '#E7298A')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF C") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfc_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#E7298A', size = 6, opacity = 0.65))

p_wrf_c
save(p_wrf_c, file = "./site_objects/wrf_c_year2.rda")

keeping in case

#save(wrfa_both, file = "./plotly_objs/wrfa_both.rda")
#save(wrfb_both, file = "./plotly_objs/wrfb_both.rda")
#save(wrfc_both, file = "./plotly_objs/wrfc_both.rda")
#save(date_vec_botha, file = "./plotly_objs/date_vec_botha.rda")
#save(date_vec_bothb, file = "./plotly_objs/date_vec_bothb.rda")
#save(date_vec_bothc, file = "./plotly_objs/date_vec_bothc.rda")
#save(both_ymina, file = "./plotly_objs/both_ymina.rda")
#save(both_ymaxa, file = "./plotly_objs/both_ymaxa.rda")

#save(both_yminb, file = "./plotly_objs/both_yminb.rda")
#save(both_ymaxb, file = "./plotly_objs/both_ymaxb.rda")

#save(both_yminc, file = "./plotly_objs/both_yminc.rda")
#save(both_ymaxc, file = "./plotly_objs/both_ymaxc.rda")